Basic concept of programming, data types, vectors and factors
Matrices and data.frames, data import/export, graphics and visualization
Tidy data, data wrangling
Today: Spatial data in R
Basic concept of programming, data types, vectors and factors
Matrices and data.frames, data import/export, graphics and visualization
Tidy data, data wrangling
Today: Spatial data in R
library(raster)
ras <- raster(ncol = 100, nrow = 100, xmn = 0, xmx = 100, ymn = 0, ymx = 100) ras
## class : RasterLayer ## dimensions : 100, 100, 10000 (nrow, ncol, ncell) ## resolution : 1, 1 (x, y) ## extent : 0, 100, 0, 100 (xmin, xmax, ymin, ymax) ## coord. ref. : NA
m <- rnorm(n = ncell(ras), mean = 0, sd = 1) values(ras) <- m plot(ras)
ras2 <- ras + 10
ras3 <- ras * ras2
ras4 <- ras^2
ras5 <- (ras - cellStats(ras, stat = mean)) /
cellStats(ras, stat = sd)
ras.na <- ras ras.na[ras.na > 0.5] <- NA
k <- matrix(c( 0, 0.5, 1,
0.5, 1.5, 2,
1.5, 6.0, 3),
ncol = 3, byrow = TRUE)
k
## [,1] [,2] [,3] ## [1,] 0.0 0.5 1 ## [2,] 0.5 1.5 2 ## [3,] 1.5 6.0 3
ras.reclass <- reclassify(x = ras, rcl = k)
ras[10, 20]
## ## 0.7534085
ras[10, 20] <- 10
ras.stack <- stack(ras, ras2) ras.stack
## class : RasterStack ## dimensions : 100, 100, 10000, 2 (nrow, ncol, ncell, nlayers) ## resolution : 1, 1 (x, y) ## extent : 0, 100, 0, 100 (xmin, xmax, ymin, ymax) ## coord. ref. : NA ## names : layer.1, layer.2 ## min values : -3.592167, 6.407833 ## max values : 10.00000, 13.60653
layer1 <- subset(x = ras.stack, subset = 1) layer1
## class : RasterLayer ## dimensions : 100, 100, 10000 (nrow, ncol, ncell) ## resolution : 1, 1 (x, y) ## extent : 0, 100, 0, 100 (xmin, xmax, ymin, ymax) ## coord. ref. : NA ## data source : in memory ## names : layer.1 ## values : -3.592167, 10 (min, max)
data <- as.data.frame(ras.stack) head(data)
## layer.1 layer.2 ## 1 -0.531517486 9.468483 ## 2 -0.183922864 9.816077 ## 3 0.003593095 10.003593 ## 4 -0.989860242 9.010140 ## 5 -3.592166533 6.407833 ## 6 -1.795983887 8.204016
elevation <- raster("Data/elevation.envi")
elevation
## class : RasterLayer ## dimensions : 709, 1306, 925954 (nrow, ncol, ncell) ## resolution : 0.008333333, 0.008333333 (x, y) ## extent : 16.90833, 27.79167, 44.33333, 50.24167 (xmin, xmax, ymin, ymax) ## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 ## data source : /Users/corneliussenf/Dropbox/Teching/Quantitative Methods/Week-13/Data/elevation.envi ## names : elevation
projection(elevation)
## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
plot(elevation)
# use stack() or brick() for multi-layer raster files
climate <- stack("Data/climate.envi")
# get the number of raster layers
nlayers(climate)
## [1] 4
# get the raster layer names names(climate)
## [1] "X.Tmax." "X.Tmin." "X.Tmean." "X.Prec."
# rename the raster layers
names(climate) <- c("Tmax", "Tmin", "Tmean", "Prec")
writeRaster(elevation,
filename= "Data/elev.envi",
overwrite = TRUE)
library(sp)
coords_srs <- data.frame(x = runif(100, 17, 27),
y = runif(100, 45, 50))
pts_srs <- SpatialPoints(coords = coords_srs)
pts_srs
## class : SpatialPoints ## features : 100 ## extent : 17.04253, 26.99759, 45.01899, 49.98391 (xmin, xmax, ymin, ymax) ## coord. ref. : NA
proj <- projection(elevation) CRS(proj)
## CRS arguments: ## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
pts_srs <- SpatialPoints(coords = coords_srs,
proj4string = CRS(proj))
plot(elevation) plot(pts_srs, add = TRUE)
library(mapview) mapview(climate)
df.climate <- extract(climate, pts_srs, df = TRUE) class(df.climate)
## [1] "data.frame"
head(df.climate)
## ID Tmax Tmin Tmean Prec ## 1 1 161.3333 60.66667 110.66666 47.83333 ## 2 2 137.1667 38.58333 87.58334 47.33333 ## 3 3 118.1667 38.00000 77.83334 56.33333 ## 4 4 154.0000 53.91667 103.58334 46.75000 ## 5 5 116.8333 28.33333 72.41666 65.33334 ## 6 6 147.9167 40.50000 93.91666 48.50000
spdf.climate <- extract(climate, pts_srs, sp = TRUE) spdf.climate
## class : SpatialPointsDataFrame ## features : 100 ## extent : 17.04253, 26.99759, 45.01899, 49.98391 (xmin, xmax, ymin, ymax) ## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 ## variables : 4 ## names : Tmax, Tmin, Tmean, Prec ## min values : 83.75, -2.5, 40.3333320617676, 43.4166679382324 ## max values : 164.16667175293, 65.6666641235352, 114.25, 86.4166641235352
You can acess SpatialPointsDataFrame the same way as data.frame!
head(spdf.climate)
## Tmax Tmin Tmean Prec ## 1 161.3333 60.66667 110.66666 47.83333 ## 2 137.1667 38.58333 87.58334 47.33333 ## 3 118.1667 38.00000 77.83334 56.33333 ## 4 154.0000 53.91667 103.58334 46.75000 ## 5 116.8333 28.33333 72.41666 65.33334 ## 6 147.9167 40.50000 93.91666 48.50000
spdf.climate$Tmax[1]
## [1] 161.3333
library(rgdal)
pts.sample <- readOGR('Data', 'pts_sample')
## OGR data source with driver: ESRI Shapefile ## Source: "Data", layer: "pts_sample" ## with 142 features ## It has 1 fields ## Integer64 fields read as strings: id
pts.sample
## class : SpatialPointsDataFrame ## features : 142 ## extent : 17.80209, 27.04407, 44.64419, 49.84834 (xmin, xmax, ymin, ymax) ## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 ## variables : 1 ## names : id ## min values : 0 ## max values : 99
writeOGR(spdf.climate, "Data", "pts_climate", driver = "ESRI Shapefile", overwrite_layer = TRUE)